Raport - badania krwi pacjentów zakazonych COVID-19

Executive summary

TODO

Użyte biblioteki

library(readxl)
library(httr)
library(arules)
library(tidyr)
library(dplyr)
library(lubridate)
library(stringr)
library(ggplot2)
library(purrr)
library(broom)
library(tibble)
library(plotly)
library(gridExtra)
library(caret)
library(survival)

Przetwarzanie zbioru danych

Odczyt danych

Odczyt danych przeprowdzony jest bezośrednio ze strony, aby zapewenić stałą wartość początkową dla każdego uruchomienia. Dodatkowo uzupełniane są brakujące dane identyfikatory pacjentów PATIENT_ID. Wartości brakujące wynikają ze sposobu przetwarzania scalonych komórek pliku źródłowego.

# https://stackoverflow.com/a/41368947
url <- "http://www.cs.put.poznan.pl/dbrzezinski/teaching/zed/wuhan_blood_sample_data_Jan_Feb_2020.xlsx"
GET(url, write_disk(tf <- tempfile(fileext = ".xlsx")))

df <- read_excel(tf) %>% fill(PATIENT_ID)

Atrybuty zbioru danych

Surowy zbiór danych zawiera 81 atrybutów i 6120 wierszy.

Atrybuty zbioru danych można potwierdzić na dwie grupy:

• Atrybuty ogólne - zawierają etykiety, podstawowe dane pacjenta(płeć, wiek, identyfikator, datę przyjęcia i wypisu, rezultat) oraz datę wykonania badania
• Wyniki badań - w kolumnach są umieszczone wartości poszczególnych badanych właściwości krwi

Brakujące dane

Poza danymi niedostępnymi (NA) w zbiorze danych wsytępują dane brakujące, oznakowane wartością -1. Poniższa tabela pokazuje ilość takich wartości w każdej kolumnie, o ile te wartości występują.

name	n
2019-nCoV nucleic acid detection	501
Platelet count	4

Dla kolumny Platelet count warkości brakujących jest mało, liczby -1 zostaną zamienione na NA.

df <- df %>% mutate(`Platelet count` = if(is.na(`Platelet count`) || `Platelet count` != -1) `Platelet count` else NA)

Dla kolumny 2019-nCoV nucleic acid detection sprawdzimy całkowitą ilość wartości różnych od NA.

name	n
2019-nCoV nucleic acid detection	501

Wszystkie wartości w kolumnie 2019-nCoV nucleic acid detection różne od NA wynoszą -1, w związku z tym ta kolumna jest właściwie nieprzydatna. Mimo to na ten moment postanowiono nie zmieniać wartości w tej kolumnie, ewentualne zmiany można wykonać w dalszych krokach.

Atrybuty ogólne

Pierwszych siedem kolumn to atrybuty ogólne, reprezentują następująco:

• PATEINT_ID - sztuczny identyfikator pacjenta
• RE_DATE - najpewniej jest to data wprowadzenia wyniku badań próbki krwi do systemu
• age - wiek pacjenta
• gender - płeć pacjenta, przyjmuje wartość 1 albo 2, ale z poisu zbioru danych nie wynika wprost, które oznaczenie jest przypisane do której płci
• Admission time - czas przyjęcia pacjenta do szpitala
• Discharge time - czas zakończenia leczenia (w skutek wypisu pacjenta lub jego śmierci)
• outcome - rezultat, przyjmuje wartość 0 jeżeli pacjent został wyleczony, lub wartość 1 jeżeli zmarł.

Nazwy wyżej wymienionych atrybutów zostaną przeformatowane, będą pisane wielkimi literami, a ewentualne spacje będą zamieniane na znak podkreślinka _, w ten sposób atrbuty ogólne będą łatwo odróżnialne od pozostałych atrybutów.

df <- df %>% rename_with(~ toupper(gsub(" ", "_", .x, fixed = TRUE)), 1:7)

Pola GENDER i OUTCOME zostaną dodatkowo zmienione na typ factor - z typu numerycznego zamienione zostaną na zmienne nominalne. Do tego trzeba przeprowadzić identyfikację typu płci.

Wartości płci nie są jawnie opisane, w artykule źródłowym została natomiast podana proporcja pacjentów danej płci, co powinno umożliwić okreslenie właściwego przypisania.

GENDER	count
1	224
2	151

Wniosek:

• Mężczyźni - wartość 1
  
• Kobiety - wartość 2

df <- df %>% 
  mutate(across("GENDER", ~factor(., levels=c(1,2), labels=c("MALE", "FEMALE")))) %>%
  mutate(across("OUTCOME", ~factor(., levels=c(0,1), labels=c("CURED","DECEASED"))))

Poniżej można zapoznać się z analizą wartości dla atrybutów ogólnych. Wartość atrybutu RE_DATE ma unikalną wartość dla kazdego wiersza.

	RE_DATE
	Min. :2020-01-10 19:45:00
	1st Qu.:2020-02-04 13:44:00
	Median :2020-02-09 12:42:30
	Mean :2020-02-08 07:00:02
	3rd Qu.:2020-02-13 10:34:00
	Max. :2020-02-18 17:49:00
	NA’s :14

Pozostałe z atrybtów ogólnych, są wspólne dla każdego pacjenta, zastosowano transformację do zniwelowania skutków duplikacji.

	PATIENT_ID	AGE	GENDER	ADMISSION_TIME	DISCHARGE_TIME	OUTCOME
	Min. : 1.0	Min. :18.00	MALE :224	Min. :2020-01-10 15:52:20	Min. :2020-01-23 09:09:23	CURED :201
	1st Qu.: 94.5	1st Qu.:46.00	FEMALE:151	1st Qu.:2020-02-01 19:27:40	1st Qu.:2020-02-11 13:39:21	DECEASED:174
	Median :188.0	Median :62.00		Median :2020-02-04 22:30:34	Median :2020-02-16 17:40:07
	Mean :188.0	Mean :58.83		Mean :2020-02-04 20:13:51	Mean :2020-02-15 16:42:59
	3rd Qu.:281.5	3rd Qu.:70.00		3rd Qu.:2020-02-10 04:11:10	3rd Qu.:2020-02-19 11:47:14
	Max. :375.0	Max. :95.00		Max. :2020-02-17 21:30:07	Max. :2020-03-04 16:21:51

Liczba unikatowych identyfikatorów pacjentów wynosi 375.

Wyniki badań

Pozostałe 74 kolumn zawiera wartości atrybutów, które reprezentują rodzaje miar zebranych podczas każdego badania.

Nazwy kolumn zostały poddane następującym transformacjom:

• Zamieniono długie nazwy kolumn na krótkie, skrótowe odpowiedniki, na podstawie nazw medycznych (częściowy code book dostępny jest w repozytorium)
• oznaczenia (%) zamieniono na suffix _percent
• oznaczenia (#) usunięto
• nazwy zaczynajace się wielką literą i poza tym skłądające się z małych liter zamieniono na odpowiedniki skłądające się wyłącznie z małych znaków, jednak skrótowce zostały zachwowane
• znak spacji i myślnika - zamieniono na podkreślenie dolne _ ze względu na łatwość odwołania się do kolumn

Poniższy blok kodu szczegółowo ukazuje zatosowane transfomrmacje, jednocześnie prezentowana są ostatenczne nazwy kolumn wyników badań.

df <- df %>%
  rename("NT-proBNP" = "Amino-terminal brain natriuretic peptide precursor(NT-proBNP)") %>%
  rename("2019-nCov-detection" = "2019-nCoV nucleic acid detection") %>%
  rename("aPPT" = "Activation of partial thromboplastin time") %>%
  rename("ALP" = "Alkaline phosphatase") %>%
  rename("AAT" = "aspartate aminotransferase") %>%
  rename("FDPs" = "Fibrin degradation products") %>%
  rename("ALT" = "glutamic-pyruvic transaminase") %>%
  rename("hs-CRP" = "High sensitivity C-reactive protein") %>%
  rename("hs-cTnI" = "Hypersensitive cardiac troponinI") %>%
  rename("HCV-abs-quant" = "HCV antibody quantification") %>%
  rename("HIV-abs-quant" = "HIV antibody quantification") %>%
  rename("INR" = "International standard ratio") %>%
  rename("LDH" = "Lactate dehydrogenase") %>%
  rename("MCH" = "mean corpuscular hemoglobin") %>%
  rename("MCHC" = "mean corpuscular hemoglobin concentration") %>%
  rename("MCV" = "mean corpuscular volume") %>%
  rename("MPV" = "Mean platelet volume") %>%
  rename("P-LCR" = "platelet large cell ratio") %>%
  rename("PDW" = "PLT distribution width") %>%
  rename("q-t-pallidum-abs" = "Quantification of Treponema pallidum antibodies") %>%
  rename("RCDW-SD" = "RBC distribution width SD") %>%
  rename("RBC_count" = "Red blood cell count") %>%
  rename("RCDW" = "Red blood cell distribution width") %>%
  rename("TNF-alfa" = "Tumor necrosis factorα") %>%
  rename("WBC_count" = "White blood cell count") %>%
  rename("gamma-GT" = "γ-glutamyl transpeptidase") %>%
  rename_with(~str_c(str_replace(.x, "\\(%\\)", ""), "_percent"), contains("(%)")) %>%
  rename_with(~str_replace(.x, "\\(#\\)", "")) %>%
  rename_with(tolower, matches("^[[:upper:]]{1}[[:lower:]]+", ignore.case = FALSE)) %>%
  rename_with(~str_replace_all(.x, "[ |-]", "_"))
str_sort(colnames(df[,-(1:7)]))

##  [1] "2019_nCov_detection"    "AAT"                    "albumin"               
##  [4] "ALP"                    "ALT"                    "antithrombin"          
##  [7] "aPPT"                   "basophil_count"         "basophil_percent"      
## [10] "calcium"                "corrected_calcium"      "creatinine"            
## [13] "D_D_dimer"              "direct_bilirubin"       "eGFR"                  
## [16] "eosinophil_count"       "eosinophils_percent"    "ESR"                   
## [19] "FDPs"                   "ferritin"               "fibrinogen"            
## [22] "gamma_GT"               "globulin"               "glucose"               
## [25] "HBsAg"                  "HCO3_"                  "HCV_abs_quant"         
## [28] "hematocrit"             "hemoglobin"             "HIV_abs_quant"         
## [31] "hs_CRP"                 "hs_cTnI"                "indirect_bilirubin"    
## [34] "INR"                    "interleukin_10"         "interleukin_1β"        
## [37] "interleukin_2_receptor" "interleukin_6"          "interleukin_8"         
## [40] "LDH"                    "lymphocyte_count"       "lymphocyte_percent"    
## [43] "MCH"                    "MCHC"                   "MCV"                   
## [46] "monocytes_count"        "monocytes_percent"      "MPV"                   
## [49] "neutrophils_count"      "neutrophils_percent"    "NT_proBNP"             
## [52] "P_LCR"                  "PDW"                    "PH_value"              
## [55] "platelet_count"         "procalcitonin"          "prothrombin_activity"  
## [58] "prothrombin_time"       "q_t_pallidum_abs"       "RBC_count"             
## [61] "RCDW"                   "RCDW_SD"                "serum_chloride"        
## [64] "serum_potassium"        "serum_sodium"           "thrombin_time"         
## [67] "thrombocytocrit"        "TNF_alfa"               "total_bilirubin"       
## [70] "total_cholesterol"      "total_protein"          "urea"                  
## [73] "uric_acid"              "WBC_count"

Czyszczenie zbioru danych

Niniejsza sekcja jest poświęcona wstępnemu przetwarzniu zbiorów danych, ze szczególnym uwzględnieniem usuwania i scalania wybranych rekordów.

Puste rekordy w zbiorze

W powyższej głównej zauważono wystąpienia wartości brakujących dla kolumny RE_DATE. Szybkie sprawdzenie wykazuje, że dla wierszy, gdzie nie ma atrybut RE_DATE przyjmuje wartość ‘NA’, wartości wszystkich atrybutów pochodzących z badań krwi, również przyjmują wartość ‘NA’.

all(sapply(df[is.na(df$RE_DATE), -c(1,3:7)], is.na))

## [1] TRUE

Fakt, że wiersze te są właściwie puste, jest podstawą do ich usunięcia ze zbioru.

df <- df[!is.na(df$RE_DATE),]

Po tej zmianie w zbiorze pozostało 6106 wierszy i 361 unikalnych identyfikatorów pacjentów.

Kategorie badań

Wstępnę analiza danych sugeruje wskazuje, że wiele wartości atrybutów jest oznaczona wartością ‘NA’. Poniższy wykresy wskazuje, ilość wystąpień rekordów badań z danę liczbą wartości (różnych od NA).

Z powyższego histogramu można wyciągnąć kilka wniosków

• Znacząca większość rekordów ma zaledwie jeden wynik badania (jedną wartość).
• Interesujący jest fakt, że druga i trzecia pod względem liczości grupa rekordów ma odpowiednio po 24 i 23 elementy, taki fakt może sugerować, że elementy mogą występować w grupach. Możliwa jest dalsza analiza współwystępowania ze sobą atrybutów, np. przy pomocy analizy zbiorów częstych.
• Ze względu na to, że właściwie nie ma rekordów, które opisują wszystkie 74 atrybuty, konieczne będzie scalanie rekordów w dalszych krokach.

Współwystępowanie kategorii

W poprzedniej sekcji stwierdzono, że można przeanalizować współwystępowanie ze sobą poszczególnych atrybutów w rekordach. W związku z tym tabela wyników badań zostanie przekształcona do zbioru binarnych transakcji.

transactions <- sapply(df[,-(1:7)], function(x) !is.na(x))

Sprawdzona zostanie ilość unikalnych transakcji.

length(unique(apply(transactions, 1, function(x) which(x))))

## [1] 176

W zbiorze 6106 transakcji jest zaledwie 176 różnych typów.

Aby sprawdzić, czy atrybuty wsytępują w grupach, wyszukiwane są maksymalne domknięte zbiory częste, dla mimalnej wartości support równej 5%.

itemsets <- apriori(transactions, parameter = list(target="closed frequent itemsets", support=.05, minlen=1, maxlen=30, maxtime=60))

## Apriori
## 
## Parameter specification:
##  confidence minval smax arem  aval originalSupport maxtime support minlen
##          NA    0.1    1 none FALSE            TRUE      60    0.05      1
##  maxlen                   target  ext
##      30 closed frequent itemsets TRUE
## 
## Algorithmic control:
##  filter tree heap memopt load sort verbose
##     0.1 TRUE TRUE  FALSE TRUE    2    TRUE
## 
## Absolute minimum support count: 305 
## 
## set item appearances ...[0 item(s)] done [0.00s].
## set transactions ...[74 item(s), 6106 transaction(s)] done [0.00s].
## sorting and recoding items ... [63 item(s)] done [0.00s].
## creating transaction tree ... done [0.00s].
## checking subsets of size 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 done [28.18s].
## filtering closed item sets ... done [19.86s].
## sorting transactions ... done [0.00s].
## writing ... [130 set(s)] done [0.35s].
## creating S4 object  ... done [1.10s].

mc_itemsets <- itemsets[is.maximal(itemsets)]

inspect(sort(mc_itemsets, by = 'support'))

##      items                     support transIdenticalToItemsets count
## [1]  {hemoglobin,                                                    
##       eosinophils_percent,                                           
##       basophil_percent,                                              
##       platelet_count,                                                
##       monocytes_percent,                                             
##       RCDW,                                                          
##       neutrophils_percent,                                           
##       MCV,                                                           
##       hematocrit,                                                    
##       WBC_count,                                                     
##       MCHC,                                                          
##       lymphocyte_count,                                              
##       RBC_count,                                                     
##       eosinophil_count,                                              
##       neutrophils_count,                                             
##       MPV,                                                           
##       RCDW_SD,                                                       
##       lymphocyte_percent,                                            
##       P_LCR,                                                         
##       monocytes_count,                                               
##       PDW,                                                           
##       basophil_count,                                                
##       MCH,                                                           
##       thrombocytocrit}      0.13609564                0.1354406   831
## [2]  {glucose}              0.12692434                0.0000000   775
## [3]  {serum_chloride,                                                
##       ALP,                                                           
##       albumin,                                                       
##       total_bilirubin,                                               
##       indirect_bilirubin,                                            
##       total_protein,                                                 
##       urea,                                                          
##       corrected_calcium,                                             
##       serum_potassium,                                               
##       direct_bilirubin,                                              
##       total_cholesterol,                                             
##       AAT,                                                           
##       uric_acid,                                                     
##       HCO3_,                                                         
##       calcium,                                                       
##       LDH,                                                           
##       globulin,                                                      
##       gamma_GT,                                                      
##       hs_CRP,                                                        
##       serum_sodium,                                                  
##       ALT,                                                           
##       eGFR,                                                          
##       creatinine}           0.11021946                0.1094006   673
## [4]  {hs_cTnI}              0.08303308                0.0000000   507
## [5]  {2019_nCov_detection}  0.08205044                0.0000000   501
## [6]  {NT_proBNP}            0.07779234                0.0000000   475
## [7]  {procalcitonin}        0.07517196                0.0000000   459
## [8]  {PH_value}             0.06288896                0.0000000   384
## [9]  {ESR}                  0.06272519                0.0000000   383
## [10] {prothrombin_time,                                              
##       antithrombin,                                                  
##       prothrombin_activity,                                          
##       fibrinogen,                                                    
##       thrombin_time,                                                 
##       D_D_dimer,                                                     
##       FDPs,                                                          
##       INR,                                                           
##       aPPT}                 0.05371765                0.0000000   328

Otrzymaliśmy 10 maksymalnych domkniętych zbiorów częstych, z czego 7 to zbiory jednolementowe. Pozostałe 3 zbiory (1, 3 i 10) zawierają wiele elementów. Na podstawie elementów składowych, zbiory wieolemelemnetowe zostały zakwalifikowane następująco:

• Zbiór 1 - zawiera elementy typowe dla badania morfologii krwi (ang. Compete blood count (CBC))
• Zbiór 3 - zawiera elementy typowe dla badania biochemii krwi (ang. Blood biochemistry)
• Zbiór 10 - zawiera elementy typowe dla badania krzepliwości krwi (ang. Coagulation)

Na podstawie uzyskanego wyniku można wnioskować, że ilość elementów w rekordzie może wynikać z procesu technologicznego analizowania próbek krwi. Sugeruje to, że możliwe będzie scalanie rekordów, np na podstawie dnia wykonania badania (jeśli jednemu pacjentowi przeprowadzono wiele badań w ciagu dnia).

Badania w czasie

W tej sekcji przeprowadzona zostanie analiza występowania badań w czasie. Ponadto rekordy zostaną pogrupowane według przynależności do zbiorów wyznaczonych w poprzedniej sekcji oraz dwóch dodatkowych klas NONE jeśli rekord nie należy do żadnej kategorii i MULTIPLE jeśli rekod należy do więcej niż jednej kategorii.

Poniższe blok kodu generuje oetykietowany zbiór rekordów badań, ograniczony jedynie do podstawowych informacji.

class_labels <- as(items(sort(mc_itemsets, by="support")), "list")
class_labels <- sapply(class_labels, function(x) x[1])
class_labels[1] <- "Complete blood count"
class_labels[3] <- "Blood biochemistry"
class_labels[10] <- "Coagulation"

classes <- as(items(sort(mc_itemsets, by="support")), "matrix")
labeled <- df %>% rowwise() %>%
  mutate(CLASS = paste(
    which(
      apply(
        (matrix(
          !is.na(c_across(-(1:7))), 
          nrow=nrow(classes), 
          ncol=ncol(classes), 
          byrow=TRUE, 
          dimnames=dimnames(classes)
          ) & classes) == classes, 
        1, 
        all)
      ), collapse=""), 
    .after=OUTCOME ) %>%
  mutate(CLASS = if(CLASS != "" && as.integer(CLASS)>10) "MULTIPLE" else CLASS) %>%
  mutate(CLASS = if(CLASS != "" && CLASS != "MULTIPLE") class_labels[as.integer(CLASS)] else CLASS) %>%
  mutate(CLASS = str_replace(CLASS, "^$", "NONE")) %>%
  mutate(CLASS = str_trim(CLASS)) %>%
  mutate(CLASS = factor(CLASS)) %>%
  rowwise() %>%
  mutate(SIZE = sum(!is.na(c_across(-c(1:8)))), .after=CLASS) %>%
  select(1:9) %>%
  ungroup()

knitr::kable(
  labeled %>% group_by(CLASS) %>% summarize(count = n(), .groups="drop") %>% arrange(desc(count))
)

CLASS	count
NONE	1433
Complete blood count	816
glucose	584
Blood biochemistry	507
2019_nCov_detection	500
MULTIPLE	485
hs_cTnI	386
PH_value	375
ESR	368
Coagulation	300
NT_proBNP	190
procalcitonin	162

Powyższa tabela prezentuje ilość rekordów przydzielonych do każdej z klas.

Poniższy wykres prezentuje badania danej kategorii według daty wykonania. Dodatkowo każdy punkt jest opisany identyfikatorem pacjenta.

Z powyższego wykresu można wywnioskować kilka rzeczy:

• Większość badań została przeprowadzona po 27 stycznia 2020, wszystkie badania zostały wykonane do 19 lutego 2020. Można rozważyć przeprowadzenie dodatkowej analizy skupiającej się na ilości przyjęć i wypisów ze szpitala, w zależności od dnia, aby wyciągnąć dodatkowe wnioski. Fakt ten moze wynikać ze zwiększonej liczby przyjmowanych pacjentów po 27 stycznia 2020.
• Część badań nie była w praktyce przeprowdzana (lub na małą skalę) przed 27 stycznia 2020
• Na wykresie można zaobserwować, że poszczególne punkty są często zgrupowane w “paski” - takie grupy rekordów cechują się małymi różnicami w wartości atrybutu RE_DATE. Biorąc pod uwagę fakt, że takie punkty najcześciej przynależą do różnych pacjentów i do tej samej kategorii, należy wnioskować, że poszczególne rodzaje badań były wykonywane równolegle dla wielu pacjentów.

Jednocześnie nie znaleziono żadnych czynników, które podważają sens scalania rekordów.

Badania pacjentów

Kolejna analiza dotyczy ilości różnych atrybutów zbadanych u każdego pacjenta. W tym celu wyznaczamy wektor wartości ostatnich badań każdego z pacjentów, tj wynik każdego badania jest najmłodszy i różny od NA, jeśli dany pacjent nie ma przypisanego żadnego rekordu zawierajacego wartość danego atrybutu, to wartość ta nadal będzie oznakowana jako NA.

Powyższy histogram obrazuje rozkład ilości badań wykonanych pacjentom. Zauważalna jest grupa wartości odstających - pacjentów którym wykonano znacząco mniej rodzajów badań. Rekordy tych pacjentów zostaną usunięte ze zbioru danych.

Usuwanie pacjentów z znacząco mniejszą liczbą badań

W celu usunięcia pacjentów z mniejszą ilością badań, zostanie wyznaczony próg minimlanej ilości badań. Próg ten przyjmuje wartość między 0 a 1, należy go interpretować jako minimalną procentową ilość badań, które musi mieć wykonany pacjent. Maksimum stanowi ilość wszystkich badanych atrybutów, czyli 74.

Wartość progu została wyznaczona analitycznie, w następujący sposób; dla każdej wartości progu z zakresu od 0.5 do 1 z krokiem 0.025, wyznaczona został procent pozostałych (nie odrzuconych) pacjentów, w stosunku do liczby początkowej, a także procent kolumn, które nie zawierają żadnej wartości NA.

Uzyskane dane zostały przedstawione na poniższym wykresie.

Wybrana wartość została oznaczona na wykresie czarną przerywaną linią, wynosi 0.65 i odpowiada 48 atrybutom. Taka wartość progowa pozwala na zachowanie znaczącej większości pacjentów w zbiorze, uzyskanie lepszego efektu poprzez jej zwiększenie, wiązałoby się ze znaczącym zmniejszeniem ilości pacjentów w zbiorze wynikowym.

patients_to_be_removed <- patient_last_results %>% 
  filter(number_of_tests < (ncol(df)-7) * threshold) %>% 
  select(PATIENT_ID)

df <- df %>% filter(PATIENT_ID %!in% patients_to_be_removed$PATIENT_ID)

Po usunięciu rekordów pacjentów, u których przeprowadzono mało badań, pozostało 354 unikalnych pacjentów.

Brakujące wartości w kategoriach

Poniższy histogram prezentuje ilość kolumn, według ilości brakujących wartości badań w każdej kolumnie, biorąc pod uwagę wszystkie wyniki każdego z pacjentów. Dodatkowo nazwy kolumn zostały oznaczone według przynależności do jednej z grup badań (sekcja “Współwystępowanie badań”), lub opsiane wartością NONE, jesli nie należą do żadnej z nich.

Na powyższym wykresie możemy zaobserwować, że znacząca część badań została wykonana wszystkim pozostałym w zbiorze pacjentom. Dodatkowo, wszystkie testy z grupy biochemii krwi i morfologii zostały przeprowadzone u większości pacjentów.

W kolejnych krokach można rozważyć odrzucenie pacjentów, którzy nie mają wykonanych wszystkich testów w tych dwóch grupach.

Agregacja wierszy według daty badania

W zwiazku z brakiem wykrytycyh preciwskazań, w tym kroku przeprowadzono scalanie wyników badań każdego pacjenta, według dnia jego wykonania, tj. wszystkie rekordy z danego dnia zostaną scalone w jeden rekord, który zawiera wszystkie, najbardziej aktualne wyniki z danego dnia.

df <- df %>% 
  mutate(RE_DATE = floor_date(RE_DATE, '1 day')) %>% 
  group_by(PATIENT_ID, RE_DATE) %>% 
  fill(everything()) %>% 
  summarise_all(last) %>% 
  ungroup()

Po tej transformacji, w zbiorze danych pozostało 1688 rekordów.

Powyższy histogram prezentuje ilość rekordów zawierających daną ilość wyników badań. Pomimo znacznej ilości rekordów zawierajacych tylko jeden wynik badania, obserwujemy wzrost w ilości rekordów zawierających ponad 40 wyników badań, w porównaniu do wyników przed agregacją.

Zbiór danych

Ta sekcja jest poświęcona analizie danych w zbiorze.

Analiza wybranych parametrów podstawowych

Płeć

Poniższy wykres prezentuje ilość przypadków z dany rezultatem leczenia w stosunku do płci pacjenta. Obserwujemy relatywnie małą ilość zmarłych kobiet w stosunku do zmarłych mężczyzn.

Wiek

Poniższy wykres przedstawia rozkład wieku w zależności od wyniku leczenia i płci pcajenta. Większość ze zmarłych pacjentów ma powyżej 60 lat, co sugeruje że grupa starszych pacjentów jest bardziej narażona na śmierć w wyniku choroby spowodowanej zarażeniem wirusem.

Analiza wartości parametrów

Poniższy zbiór wykresów przedstawia rozkład poszczególnych wyników badań, w zależności od rezultatu leczenia. Część wykresów jest nieczytelna przez występowanie wartości odstających (outlierów).

W niektórych wypadkach możemy zaobserwować istotne różnice między zbiorami pacjentów wyleczonych i zmarłych, są to między innymi parametry takie jak: albumin, hs_CRP, lymphoocytes_percent i monocythes_percent.

Korelacja

Poniższy wykres przedstawia wartość współczynnika korelacji Pearsona między wszystkimi parametrami liczbowymi w zbiorze. Parametry nominalne GENDER i OUTCOME, ze względu na to zę mają tylko dwie wartości, zostały zakodowane jako isMALE i isCURED. Dodatkowo, ze parametr 2019_nCov_detection ze względu na stałą wartość, został zamieniony na parametr has_nCOV_test, który indykuje czy w danym badaniu wykonano taki test.

Ponadto w poniższej tabeli zebrano 30 par atrybutów z największym współczynnikiem korelacji Pearsona.

rowname	colname	value
MPV	P_LCR	0.99
INR	prothrombin_time	0.99
platelet_count	thrombocytocrit	0.98
HBsAg	interleukin_6	0.98
direct_bilirubin	total_bilirubin	0.98
lymphocyte_percent	neutrophils_percent	-0.98
procalcitonin	q_t_pallidum_abs	0.95
D_D_dimer	FDPs	0.95
P_LCR	PDW	0.95
MPV	PDW	0.94
ALT	interleukin_1β	0.93
lymphocyte_count	monocytes_count	0.88
serum_chloride	serum_sodium	0.86
eosinophil_count	eosinophils_percent	0.86
AAT	interleukin_1β	0.86
hematocrit	hemoglobin	0.84
MCH	MCV	0.82
indirect_bilirubin	total_bilirubin	0.80
interleukin_6	interleukin_8	0.80
RCDW	RCDW_SD	0.79
aPPT	prothrombin_time	0.79
monocytes_percent	neutrophils_percent	-0.77
eGFR	urea	-0.77
creatinine	urea	0.75
albumin	calcium	0.74
isCURED	neutrophils_percent	-0.73
isCURED	lymphocyte_percent	0.72
neutrophils_count	neutrophils_percent	0.71
albumin	total_protein	0.70
lymphocyte_percent	neutrophils_count	-0.70

Można zaobserwować wiele mocno skorelowanych parametrów, m.in:

• miary MPV, P-LCR i PDW odnoszą się do wielkości płytek krwi (odpowiednio: mean platelet volume, platelet - large cell ratio i platelet distribution width), więc ich wysokie skorelowanie jest łatwe do wyjaśnienia.
• lymphocytes_percent i neutrophils_percent - określają proporcje białych krwniek danego typu w krwi pacjenta
• INR z definicji jest właściwie liniowo zależna od prothrombin_time
• direct_bilirubin jest składową total_bilirubin

Dodatkowo sprawdzona zostanie korelacja poszczególnych parametrów z wynikiem leczenia chorego (isCured). Poniższa tabela zawiera wartości korelacji Pearsona wraz z przedziałem ufności i p-wartością (wszystkie z pokazanych korelacji są istotne statysytycznie) dla 15 najmocniej skorelowanych parametrów.

attribute	estimate	p.value	conf.low	conf.high
neutrophils_percent	-0.7301461	0	-0.7586392	-0.6988654
lymphocyte_percent	0.7234883	0	0.6915912	0.7525691
albumin	0.6880859	0	0.6524294	0.7207028
prothrombin_activity	0.6547754	0	0.6084767	0.6966318
hs_CRP	-0.6509896	0	-0.6910468	-0.6069460
D_D_dimer	-0.6376192	0	-0.6819778	-0.5885867
LDH	-0.6270327	0	-0.6648056	-0.5860615
neutrophils_count	-0.6083675	0	-0.6470960	-0.5665074
FDPs	-0.6042409	0	-0.6689583	-0.5304310
AGE	-0.5761382	0	-0.6071595	-0.5433633
calcium	0.5703332	0	0.5257539	0.6117881
platelet_count	0.5419441	0	0.4952125	0.5855487
monocytes_percent	0.5415337	0	0.4947996	0.5851444
eGFR	0.4909936	0	0.4402769	0.5385870
urea	-0.4841428	0	-0.5321758	-0.4330031

Z powyższej tabli wynika, że wpływ na wynik choroby pacjenta mają:

• stan układu odpornościowego (wysoka proporcja lymphocyte i niska proporcja neutrophils)
• stan układu krzepnięcia (prothrombin_activity, D_D_dimer, FDPs)
• hs_CRP i LDH których podwyższony stan wskazuje na dużą ilość uszkodzonych komórek
• wysoki poziom albumin - białek odpowiadających za transport odczynników we krwi (w tym leków) i regulację ilości wody we krwi
• wiek

Analiza przyjęć i wypisów

Poniższy wykres pokazuje ilość pacjentów przyjętych i wypisanych danego dnia ze szpitala, z podziałem na wynik leczenia.

Analiza przyjęć i wypisów nie przynosi żadnych istotnych wniosków. Obserwujemy jedynie zwiększoną ilość wypisów wyleczonych pacjentów między 16 a 20 lutego. Obserwowane luki w czasach przyjęć i wypisów mogą wynikać z podziału zbioru danych na część treningową (analizowaną) i testową (wyłączoną z niniejszej analizy).

Czas hosptializacji

Poniżesz wykresy prezentują odpowiednio rozkład czasu pobytu pacjentów w szpitalu w zależności od wyników leczenia a także histogram długości pobytu w szpitalu z podziałem na wynik leczenia.

Obserwujemy, że w grupy pacjentów, którzy zmarli, czas pobytu w szpitalu wynosił do 10 dni (75% przypadków). Tę cechę można potencjalnie wykorzystać przy przygotowaniu klasyfikatora.

Jednocześnie poczyniona obserwacja każe się zastanowić nad dokładnością klasyfikatora zaproponowanego przez autorów artykułu, z którego pochodzą dane. Badacze twierdzą, że dokładność wynosi ponad 90% nawet 10 dni przed poznaniem dokłądnego wyniku leczenia. Z drugiej strony, odrzucenie badań które powstały później, niż 10 dni przed pozniem ostatecznego wyniku, skutkuje odrzuceniem danych o większości zmarłych pacjentów, co pozostawia niezbalansowany zbiór danych. Jak wiadomo dokładność klasyfikatora nie jest odpowiednią do zastosowania metryką przy klasyfikacji niezbalansowanych zbiorów, a wartość miary Kappa nie została podana w treści artykułu, co utrudnia ocenienie jakości zaprezentowanego w artykule modelu.

Klasyfikator

Niniejsza sekcja opisuje proces tworzenia klasyfikatora, którego zadaniem jest ocena, czy dany pacjent przeżył chorobę spowodowaną zakażeniem COVID-19.

Dla tego celu, zbiór wyników badań zostanie zredukowany do danych pacjentów, zawierając ostatni, tj. najbardziej aktualny wynik badania. Wcześniej, ze zbioru zostaną usunięte te wyniki badań, które zostały uzyskane przed wypisem pacjenta ze szpitala. Dodatkowo usunięte zostanę atrybuty PATIENT_ID, ADMISSION_TIME, DISCHARGE_TIME i RE_DATE oraz 2019_nCOV-detection (ze względu na stałe wartości).

patients <- df %>% 
  filter(RE_DATE < DISCHARGE_TIME) %>%
  group_by(PATIENT_ID) %>% 
  fill(everything()) %>% 
  summarise_all(last) %>% 
  ungroup() %>%
  select(-PATIENT_ID, -ADMISSION_TIME, -DISCHARGE_TIME, -RE_DATE) %>%
  select(-'2019_nCov_detection')

Wybór cech

W związku z tym, że dla części pacjentów nie ma dostępnych wszystkich wyników badań, niezbędne jest dalsze przetwarzanie zbioru danych.

Grupa cech “Biochemia krwi” (sekcja Współwystępowanie kategorii) została określona jako minimalny zbiór cech, które powinny być obecne w każdym rekordzie. Ten wybór jest uzasadniony tym, że autorzy artykułu źródłowego wśrod najbardziej znaczących cech wskazali LDH i hs_CRP, które należą do tej grupy.

Po odrzuceniu rekordów, które nie zawieraływ szystkich wskazanych cech, usunięto kolumny, które zawierały wartości NA. W wyniku tego przetwarzania uzyskano czysty zbiór, który może być wykorzystany do trenowania klasyfikatora.

features <- class_labels[[3]]

patients <- patients %>% filter(complete.cases(patients %>% select(all_of(features))))

to.drop <- patients %>% 
  pivot_longer(4:76) %>% 
  mutate(value=is.na(value)) %>% 
  group_by(name) %>% 
  summarize(NAs = sum(value), .groups="drop") %>% 
  filter(NAs != 0)
#to.drop %>% arrange(NAs)

patients <- patients %>% select(-to.drop$name)

all(complete.cases(patients))

## [1] TRUE

Ostateczna wersja zbioru zawiera 345 rekordów pacjentów. Poniższa lista zawier natomiast nazwy wszystkich 43 cech wykorzystanych do klasyfikacji.

##  [1] "AGE"                 "GENDER"              "hemoglobin"         
##  [4] "serum_chloride"      "eosinophils_percent" "ALP"                
##  [7] "albumin"             "basophil_percent"    "total_bilirubin"    
## [10] "platelet_count"      "monocytes_percent"   "indirect_bilirubin" 
## [13] "neutrophils_percent" "total_protein"       "MCV"                
## [16] "hematocrit"          "WBC_count"           "MCHC"               
## [19] "urea"                "lymphocyte_count"    "RBC_count"          
## [22] "eosinophil_count"    "corrected_calcium"   "serum_potassium"    
## [25] "neutrophils_count"   "direct_bilirubin"    "lymphocyte_percent" 
## [28] "total_cholesterol"   "AAT"                 "uric_acid"          
## [31] "HCO3_"               "calcium"             "LDH"                
## [34] "monocytes_count"     "globulin"            "gamma_GT"           
## [37] "basophil_count"      "MCH"                 "hs_CRP"             
## [40] "serum_sodium"        "ALT"                 "eGFR"               
## [43] "creatinine"

Zbalansowanie zbioru prezentuje się następująco

##    CURED DECEASED 
##      191      154

W kolejnym kroku zbiór został podzielony na testowy i treningowy, odkładając 25% danych do testowania klasyfikatorów.

set.seed(23)
inTraining <- createDataPartition(
        y = patients$OUTCOME,
        p = .75,
        list = FALSE)

patientsTrain <- patients[ inTraining[,1],]
patientsTest  <- patients[-inTraining[,1],]

Klaysfikatory

W tej części zostaną wytrenowane i przetestowane dwa klasyfikatory - Random Forest i SVM. W trenowaniu zostanie wykorzystana powtarzalna walidacja krzyżowa, z podziałem na 10 części i trzema powtórzeniami.

Random forest

W modelu Random Forest optymlizowany będzie parametr mtree określający ilość parametrów, które trafią do drzewa w każdej rundzie podziałów. Jednocześnie ustalono, że wartość parametru określającą ilość drzew w lesie ntree będzie wynosiła 30.

set.seed(23)
tc <- trainControl(method="repeatedcv", number=10, repeats=3)
rf.grid <- expand.grid(mtry = 5:30)
rf.fit <- train(OUTCOME ~ .,
             data = patientsTrain,
             method = "rf",
             preProc = c("center", "scale"),
             trControl = tc,
             tuneGrid = rf.grid,
             ntree = 30)
rf.fit

## Random Forest 
## 
## 260 samples
##  43 predictor
##   2 classes: 'CURED', 'DECEASED' 
## 
## Pre-processing: centered (43), scaled (43) 
## Resampling: Cross-Validated (10 fold, repeated 3 times) 
## Summary of sample sizes: 234, 234, 235, 233, 234, 234, ... 
## Resampling results across tuning parameters:
## 
##   mtry  Accuracy   Kappa    
##    5    0.9783400  0.9561484
##    6    0.9796657  0.9588794
##    7    0.9797170  0.9591510
##    8    0.9809003  0.9613794
##    9    0.9860323  0.9718434
##   10    0.9834682  0.9666227
##   11    0.9847502  0.9692488
##   12    0.9834207  0.9665154
##   13    0.9859848  0.9718002
##   14    0.9795745  0.9585735
##   15    0.9835195  0.9666743
##   16    0.9834207  0.9665467
##   17    0.9859848  0.9717693
##   18    0.9847028  0.9691417
##   19    0.9808566  0.9613880
##   20    0.9834207  0.9666104
##   21    0.9848015  0.9693949
##   22    0.9821861  0.9638227
##   23    0.9808566  0.9611002
##   24    0.9795745  0.9586660
##   25    0.9796220  0.9586125
##   26    0.9783400  0.9560791
##   27    0.9808566  0.9612299
##   28    0.9782925  0.9560051
##   29    0.9770579  0.9534538
##   30    0.9783400  0.9560791
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 9.

Wyniki uzyskane w trakcie procesu trenowania charakteryzują się wysoką dokładnością i wartością miary Kappa, co suegeruje że klasyfikator moze dobrze wykonywać swoje zadanie. W kolejnych krokach te wyniki zostaną zweryfikowane na zbiorze testowym.

set.seed(23)
rf.predict <- predict(rf.fit, newdata=patientsTest)
confusionMatrix(data = rf.predict, patientsTest$OUTCOME)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction CURED DECEASED
##   CURED       45        2
##   DECEASED     2       36
##                                          
##                Accuracy : 0.9529         
##                  95% CI : (0.8839, 0.987)
##     No Information Rate : 0.5529         
##     P-Value [Acc > NIR] : <2e-16         
##                                          
##                   Kappa : 0.9048         
##                                          
##  Mcnemar's Test P-Value : 1              
##                                          
##             Sensitivity : 0.9574         
##             Specificity : 0.9474         
##          Pos Pred Value : 0.9574         
##          Neg Pred Value : 0.9474         
##              Prevalence : 0.5529         
##          Detection Rate : 0.5294         
##    Detection Prevalence : 0.5529         
##       Balanced Accuracy : 0.9524         
##                                          
##        'Positive' Class : CURED          
## 

Klasyfikator pomylił się w 4 przypadkach (na 85) osiągając dokładność na poziomie 95%.

W związkuz tym, że modeł Random Forest jest łatwo interpretowalny, w łatwy sposób można uzyskać ważność cech, która jest przedstawiona w tabeli poniżej.

	Overall
LDH	25.19929
neutrophils_percent	24.16844
lymphocyte_percent	16.89785
hs_CRP	16.28050
lymphocyte_count	10.97626

Należy zauważyć, że cechy LDH i hs_CRP, a także lymphocyte_percent są ważnymi czynnikami w wytrenowanym klasyfikatorze, co jest w zgodzie z wnioskami autorów ze źródłowego artykułu.

SVM

Dla modelu Support Vector Machine z kernelem liniowym, optymalizowany jest prametr kosztu.

set.seed(23)
# svmLinear2
svm.grid <- expand.grid(cost = c(.25, .5, .75, 1, 1.25))
svm.fit <- train(OUTCOME ~ .,
             data = patientsTrain,
             method = "svmLinear2",
             preProc = c("center", "scale"),
             trControl = tc,
             tuneGrid = svm.grid
             )
svm.fit

## Support Vector Machines with Linear Kernel 
## 
## 260 samples
##  43 predictor
##   2 classes: 'CURED', 'DECEASED' 
## 
## Pre-processing: centered (43), scaled (43) 
## Resampling: Cross-Validated (10 fold, repeated 3 times) 
## Summary of sample sizes: 234, 234, 235, 233, 234, 234, ... 
## Resampling results across tuning parameters:
## 
##   cost  Accuracy   Kappa    
##   0.25  0.9706401  0.9402500
##   0.50  0.9720722  0.9433396
##   0.75  0.9696068  0.9383617
##   1.00  0.9592953  0.9172186
##   1.25  0.9528851  0.9042441
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was cost = 0.5.

Podobnie jak model Random Forest, model SVM również osiąga wysokie wyniki w klasyfikacji na etapie treningu.

set.seed(23)
svm.predict <- predict(svm.fit, newdata=patientsTest)
confusionMatrix(data = svm.predict, patientsTest$OUTCOME)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction CURED DECEASED
##   CURED       45        1
##   DECEASED     2       37
##                                           
##                Accuracy : 0.9647          
##                  95% CI : (0.9003, 0.9927)
##     No Information Rate : 0.5529          
##     P-Value [Acc > NIR] : <2e-16          
##                                           
##                   Kappa : 0.9288          
##                                           
##  Mcnemar's Test P-Value : 1               
##                                           
##             Sensitivity : 0.9574          
##             Specificity : 0.9737          
##          Pos Pred Value : 0.9783          
##          Neg Pred Value : 0.9487          
##              Prevalence : 0.5529          
##          Detection Rate : 0.5294          
##    Detection Prevalence : 0.5412          
##       Balanced Accuracy : 0.9656          
##                                           
##        'Positive' Class : CURED           
## 

Na etapie testu model SVM popełnił jedynie 3 błędy (ponad 96% dokładności), przez co można go uznać za minimalnie lepszy od modelu Random Forest.

Analiza przeżycia

W tej sekcji wyniki badań pacjentów zostaną poddane analizie metodą regresji Coxa.

Zbiór danych został poddany następującym zmianom:

• usunięto wyniki badań, które mają RE_DATE późniejszy niż DISCHARGE_TIME
• dodano kolumnę TIME która określa liczbę dni, które upłynęły od dokonania badania, do uzyskania rezultatu leczenia (na potrzeby modelu Coxa)
• zakodowano zmienne kategoryczne OUTCOME i GENDER jako wartości numeryczne
• usunięto zbędne dla analizy kolumny PATIENT_ID, ADMISSION_TIME, DISCHARGE_TIME, RE_DATE oraz kolumnę 2019_nCov_detection ze względu na identyczne wartości

W kolejnych krokach dla każdego z atrybutów (niezależnie) wykonano analizę tzw. Hazard Ratio (HR) przy pomocy modelu proporcjonalnego hazardu Coxa. Hazard Ratio jest wyjaśnione następująco źródło:

Put another way, a hazard ratio above 1 indicates a covariate that is positively associated with the event probability, and thus negatively associated with the length of survival.

Poniższa tabela prezentuje czynniki zwiększające ryzyko śmierci pacjenta. Odfiltrowane zostały czynniki nieistotne statysytcznie (p-wartość > 0.05) oraz te, gdzie parametr HR nie przekroczył wartości 1.1.

	beta	HR	wald.test	p.value
basophil_count	15.9100	8.128e+06	47.32	0.0000000
HCV_abs_quant	0.6568	1.929e+00	6.34	0.0118137
MPV	0.5709	1.770e+00	162.75	0.0000000
serum_potassium	0.2852	1.330e+00	63.53	0.0000000
PH_value	0.2822	1.326e+00	5.36	0.0205911
INR	0.2714	1.312e+00	96.76	0.0000000
RCDW	0.2019	1.224e+00	130.65	0.0000000
PDW	0.1923	1.212e+00	145.21	0.0000000
neutrophils_count	0.1067	1.113e+00	369.89	0.0000000
neutrophils_percent	0.1022	1.108e+00	259.56	0.0000000

Parametrem, który w największym stopniu zwiększa ryzyko śmierci, według przygotowanego modelu jest basophil_count. Jednocześnie uzskany wynik zdaje się być wątpliwy, gdyż basophil_count ma wartość parametru HR o sześć rzędów wielkości większą, niż kolejny parametr w rankingu.

Pozostałe złe prognostyki, dotyczą między innymi zwiększonej obecności neutrofilów, krzepliwości krwi (wyskoie INR) a także jakość komórek krwi (MPV, RCDW, PDW).

Poniższa tabela przedstawia dobre prognostyki jeśli chodzi o skutki leczenia. Ponownie odlitrowano czynniki nieistotne statystycznie (p-wartość > 0.05), oraz pozostawion wartości, gdzie parametr HR ma wartość poniżej 0.9.

	beta	HR	wald.test	p.value
HCO3_	-0.1079	0.8977000	107.45	0e+00
albumin	-0.1363	0.8726000	288.36	0e+00
lymphocyte_percent	-0.1455	0.8646000	250.00	0e+00
monocytes_percent	-0.2517	0.7775000	199.03	0e+00
total_cholesterol	-0.5328	0.5870000	83.96	0e+00
GENDER	-0.8068	0.4463000	81.22	0e+00
eosinophils_percent	-1.1360	0.3211000	71.02	0e+00
basophil_percent	-1.5800	0.2059000	28.46	1e-07
lymphocyte_count	-2.1470	0.1169000	224.40	0e+00
calcium	-3.6360	0.0263700	250.16	0e+00
eosinophil_count	-6.6720	0.0012660	26.65	2e-07
thrombocytocrit	-8.6920	0.0001678	142.03	0e+00

Do pozytywnych prognostyków można zaliczyć:

• większy poziom albumin
• większą proporcję monocytów, leukocytów i eozynofilów we krwi
• płeć - kobiety mają większą szansę przeżycia
• większy cholesterol

Wśród czynników obniżających wpółczynnik ryzyka, znajduje się też basophil_percent co wydaje się być szczególnie intrygujące w kontekście wykrytych negatywnych prognostyków.

Wśród wykrytych pozytywnych prognostyków, 3 ostatnie wydają się mieć bardzo istotny wpływ na obniżenie ryzyka śmierci pacjenta.

Techniczne

Ninejsza sekcja zawiera dodatkowe informacje o technicznych parametrach raportu

Wszytkie użyte pakiety

unique(c(.packages(), loadedNamespaces()))

##  [1] "survival"     "caret"        "lattice"      "gridExtra"    "plotly"      
##  [6] "tibble"       "broom"        "purrr"        "ggplot2"      "stringr"     
## [11] "lubridate"    "dplyr"        "tidyr"        "arules"       "Matrix"      
## [16] "httr"         "readxl"       "stats"        "graphics"     "grDevices"   
## [21] "utils"        "datasets"     "methods"      "base"         "Rcpp"        
## [26] "class"        "digest"       "ipred"        "foreach"      "R6"          
## [31] "cellranger"   "plyr"         "backports"    "stats4"       "e1071"       
## [36] "evaluate"     "highr"        "pillar"       "rlang"        "curl"        
## [41] "lazyeval"     "data.table"   "rpart"        "rmarkdown"    "labeling"    
## [46] "splines"      "gower"        "htmlwidgets"  "munsell"      "compiler"    
## [51] "xfun"         "pkgconfig"    "htmltools"    "nnet"         "tidyselect"  
## [56] "prodlim"      "codetools"    "randomForest" "viridisLite"  "crayon"      
## [61] "withr"        "MASS"         "recipes"      "ModelMetrics" "grid"        
## [66] "nlme"         "jsonlite"     "gtable"       "lifecycle"    "magrittr"    
## [71] "pROC"         "scales"       "stringi"      "farver"       "reshape2"    
## [76] "timeDate"     "ellipsis"     "generics"     "vctrs"        "lava"        
## [81] "iterators"    "tools"        "glue"         "crosstalk"    "yaml"        
## [86] "colorspace"   "knitr"